dweibull (Double Weibull) Distribution

dweibull (Double Weibull) Distribution#

This notebook is a self-contained, math-first tour of SciPy’s scipy.stats.dweibull distribution.

Goals

  • Understand the shape parameter and how it controls peakedness, bimodality, and tail behavior.

  • Derive key results (PDF/CDF, moments, likelihood) with clean substitutions.

  • Implement NumPy-only sampling and validate it against SciPy.

  • See how dweibull shows up in testing, Bayesian inference, and generative/noise models.

Throughout, we use the standardized distribution (default loc=0, scale=1) unless stated otherwise.

import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio

from scipy.special import gamma
from scipy.stats import dweibull, kstest

pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)

rng = np.random.default_rng(42)

1) Title & Classification#

  • Name (SciPy): dweibull (double Weibull distribution)

  • Type: Continuous

  • Support: \(x \in (-\infty, \infty)\)

  • Parameter space (standard form): shape \(c > 0\)

  • SciPy parameterization: dweibull(c, loc=0, scale=1) with

    • \(c > 0\) (shape)

    • \(\text{loc} \in \mathbb{R}\) (location shift)

    • \(\text{scale} > 0\) (scale stretch)

2) Intuition & Motivation#

A simple generative story#

A convenient way to think about dweibull is:

  1. Draw a magnitude \(Y \ge 0\) from a (one-sided) Weibull distribution.

  2. Flip a fair coin for the sign \(S \in \{-1, +1\}\).

  3. Set \(X = S\,Y\).

This immediately explains why the distribution is symmetric about 0 and why \(|X|\) is Weibull.

What it models#

dweibull is useful when you want a symmetric distribution whose shape can morph between:

  • sharp peak at 0 with heavy (stretched-exponential) tails (\(0 < c < 1\))

  • Laplace / double-exponential (\(c = 1\))

  • bimodal shapes with a dip at 0 (\(c > 1\))

The key surprise is the last bullet: for \(c>1\), the PDF at 0 becomes zero, and the distribution has two symmetric modes away from 0.

Real-world use cases#

  • Signed magnitudes: deviations that come with a size (Weibull-like) and a random sign (e.g., anomaly sizes, symmetric measurement deviations).

  • Flexible error/noise models: use \(c\) to tune tail heaviness vs concentration near 0 (especially for \(c\le 1\)).

  • Bimodal symmetric data: when values tend to avoid 0 but cluster around \(\pm m\) for some magnitude.

Relations to other distributions#

  • If \(X \sim \texttt{dweibull}(c)\), then \(|X|\) is Weibull with the same shape parameter \(c\).

  • \(c=1\) gives the Laplace distribution: \(f(x)=\tfrac12 e^{-|x|}\).

  • \(c=2\) implies \(|X|\) is Rayleigh (with a particular scale), so dweibull becomes a symmetric “signed Rayleigh magnitude” model.

3) Formal Definition#

PDF (standardized)#

For shape \(c>0\), the probability density function is

\[ f(x\mid c) = \frac{c}{2}\,|x|^{c-1}\,\exp\left(-|x|^c\right),\qquad x\in\mathbb{R}. \]

CDF (standardized)#

The CDF has a clean piecewise form. For \(x<0\),

\[ F(x\mid c)=\tfrac12\exp\left(-|x|^c\right), \]

and for \(x\ge 0\),

\[ F(x\mid c)=1-\tfrac12\exp\left(-x^c\right). \]

Location/scale form#

SciPy’s loc and scale apply the standard transformation

\[ X = \text{loc} + \text{scale}\cdot Z,\qquad Z\sim \texttt{dweibull}(c). \]

Then

\[ f_X(x) = \frac{1}{\text{scale}}\,f_Z\!\left(\frac{x-\text{loc}}{\text{scale}}\right). \]

Quantile function (PPF)#

Because the CDF is explicit, inverse-CDF sampling is easy. For \(q\in(0,1)\),

\[\begin{split} \operatorname{PPF}(q)=\begin{cases} \text{loc} - \text{scale}\,\big[-\ln(2q)\big]^{1/c}, & 0<q<\tfrac12,\\[4pt] \text{loc} + \text{scale}\,\big[-\ln\big(2(1-q)\big)\big]^{1/c}, & \tfrac12\le q<1. \end{cases} \end{split}\]
def _validate_params(c: float, scale: float) -> None:
    if not np.isfinite(c) or c <= 0:
        raise ValueError(f"shape c must be > 0, got {c!r}")
    if not np.isfinite(scale) or scale <= 0:
        raise ValueError(f"scale must be > 0, got {scale!r}")


def dweibull_pdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    # NumPy implementation of the PDF (with loc/scale).
    _validate_params(c, scale)
    x = np.asarray(x, dtype=float)
    z = (x - loc) / scale
    az = np.abs(z)
    return (c / (2 * scale)) * np.power(az, c - 1) * np.exp(-np.power(az, c))


def dweibull_logpdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    # NumPy implementation of log-PDF (handles x==loc explicitly).
    _validate_params(c, scale)
    x = np.asarray(x, dtype=float)
    z = (x - loc) / scale
    az = np.abs(z)

    out = np.empty_like(az)
    zero = az == 0

    # At z=0: pdf(0) is 0 if c>1, 1/(2*scale) if c=1, and +inf if c<1.
    if c > 1:
        out[zero] = -np.inf
    elif np.isclose(c, 1.0):
        out[zero] = -np.log(2 * scale)
    else:
        out[zero] = np.inf

    nz = ~zero
    out[nz] = (
        np.log(c)
        - np.log(2 * scale)
        + (c - 1) * np.log(az[nz])
        - np.power(az[nz], c)
    )
    return out


def dweibull_cdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
    # NumPy implementation of the CDF using expm1 for precision near 0.
    _validate_params(c, scale)
    x = np.asarray(x, dtype=float)
    z = (x - loc) / scale

    out = np.empty_like(z)
    neg = z < 0
    pos = ~neg

    zneg = -z[neg]
    zpos = z[pos]

    out[neg] = 0.5 * (1.0 + np.expm1(-np.power(zneg, c)))
    out[pos] = 0.5 - 0.5 * np.expm1(-np.power(zpos, c))

    return out


def dweibull_ppf(q, c: float, loc: float = 0.0, scale: float = 1.0):
    # Quantile function (inverse CDF).
    _validate_params(c, scale)
    q = np.asarray(q, dtype=float)

    out = np.empty_like(q)
    out[q <= 0] = -np.inf
    out[q >= 1] = np.inf

    mid = (q > 0) & (q < 1)
    qmid = q[mid]

    left = qmid < 0.5
    right = ~left

    out_mid = np.empty_like(qmid)
    out_mid[left] = -np.power(-np.log(2.0 * qmid[left]), 1.0 / c)
    out_mid[right] = np.power(-np.log(2.0 * (1.0 - qmid[right])), 1.0 / c)

    out[mid] = loc + scale * out_mid
    return out

# Quick sanity check: CDF(PPF(q)) ≈ q
c_test = 1.3
q_grid = np.linspace(1e-6, 1 - 1e-6, 2000)
x_from_q = dweibull_ppf(q_grid, c=c_test)
q_back = dweibull_cdf(x_from_q, c=c_test)

max_err = np.max(np.abs(q_back - q_grid))
max_err

1.1102230246251565e-16

4) Moments & Properties#

Absolute moments (key identity)#

For the standardized distribution \(Z\sim\texttt{dweibull}(c)\), a very useful identity is

\[ \mathbb{E}[|Z|^k] = \Gamma\!\left(1+\frac{k}{c}\right),\qquad k>-c. \]

In particular, all positive moments exist for any \(c>0\).

Mean, variance, skewness, kurtosis#

Because the PDF is symmetric, all odd central moments are 0 (when they exist). For \(c>0\):

  • Mean: \(\mathbb{E}[Z]=0\) and \(\mathbb{E}[X]=\text{loc}\).

  • Variance: \(\operatorname{Var}(Z)=\Gamma\!\left(1+\frac{2}{c}\right)\), so \(\operatorname{Var}(X)=\text{scale}^2\,\Gamma\!\left(1+\frac{2}{c}\right)\).

  • Skewness: 0.

  • Excess kurtosis:

\[ \gamma_2 = \frac{\Gamma\!\left(1+\frac{4}{c}\right)}{\Gamma\!\left(1+\frac{2}{c}\right)^2}-3. \]

MGF / characteristic function#

  • The characteristic function \(\varphi(t)=\mathbb{E}[e^{itZ}]\) always exists. Because of symmetry, \(\varphi(t)=\mathbb{E}[\cos(tZ)]\).

  • The MGF \(M(t)=\mathbb{E}[e^{tZ}]\) depends on \(c\):

    • \(c>1\): exists for all real \(t\) (tails decay faster than exponential).

    • \(c=1\): exists only for \(|t|<1\) (Laplace case).

    • \(0<c<1\): diverges for any \(t\ne 0\) (stretched-exponential tails).

A useful analytic representation is the even-moment series (when it converges):

\[ M(t)=\sum_{n=0}^\infty \frac{t^{2n}}{(2n)!}\,\Gamma\!\left(1+\frac{2n}{c}\right). \]

Entropy#

The differential entropy for the standardized distribution is

\[ H(Z)=1-\ln c + \ln 2 + \gamma\,\Big(1-\frac{1}{c}\Big), \]

where \(\gamma\approx 0.57721\) is the Euler–Mascheroni constant. With scaling, \(H(X)=H(Z)+\ln(\text{scale})\).

Modes#

For \(c\le 1\), the distribution is unimodal with a peak at 0 (in fact, the PDF is infinite at 0 when \(c<1\)). For \(c>1\), the PDF at 0 is 0 and there are two modes at

\[ \pm\left(\frac{c-1}{c}\right)^{1/c}. \]
def dweibull_theoretical_stats(c: float, loc: float = 0.0, scale: float = 1.0):
    # Return mean, var, skewness, excess kurtosis for dweibull(c, loc, scale).
    _validate_params(c, scale)

    mean = loc
    var = (scale**2) * gamma(1.0 + 2.0 / c)

    # Symmetry => skewness = 0
    skew = 0.0

    m2 = gamma(1.0 + 2.0 / c)
    m4 = gamma(1.0 + 4.0 / c)
    excess_kurtosis = m4 / (m2**2) - 3.0

    return mean, var, skew, excess_kurtosis


def dweibull_entropy(c: float, scale: float = 1.0):
    # Differential entropy for dweibull(c, loc=0, scale).
    _validate_params(c, scale)
    return 1.0 - np.log(c) + np.log(2.0 * scale) + np.euler_gamma * (1.0 - 1.0 / c)


c_demo = 0.7
mean_th, var_th, skew_th, kurt_th = dweibull_theoretical_stats(c_demo)
H_th = dweibull_entropy(c_demo)

mean_sp, var_sp, skew_sp, kurt_sp = dweibull.stats(c_demo, moments="mvsk")
H_sp = dweibull.entropy(c_demo)

(
    np.array([mean_th, var_th, skew_th, kurt_th]),
    np.array([mean_sp, var_sp, skew_sp, kurt_sp]),
    float(H_th),
    float(H_sp),
)

(array([ 0.      ,  5.029145,  0.      , 13.777367]),
 array([ 0.      ,  5.029145,  0.      , 13.777367]),
 1.802443982398021,
 1.802443982398021)

5) Parameter Interpretation (How Shape Changes)#

Shape parameter \(c\)#

The single shape parameter controls multiple behaviors at once:

  • Near 0: the factor \(|x|^{c-1}\) decides what happens at the origin.

    • \(c<1\): \(|x|^{c-1}\to\infty\) → infinite spike at 0.

    • \(c=1\): finite value at 0 (Laplace).

    • \(c>1\): \(|x|^{c-1}\to 0\) → density drops to 0 at 0 (bimodal).

  • Tails: \(\exp(-|x|^c)\) controls tail decay.

    • Smaller \(c\) → heavier (slower) tail decay.

    • Larger \(c\) → lighter (faster) tail decay.

loc and scale#

  • loc shifts the distribution left/right (median and mean move to loc).

  • scale stretches the distribution; variance scales like \(\text{scale}^2\).

# PDF shapes for different c
x = np.linspace(-4, 4, 2000)

c_values = [0.5, 0.8, 1.0, 1.5, 3.0]
fig = go.Figure()

for c in c_values:
    y = dweibull_pdf(x, c)

    # For c<1, the PDF spikes to +inf at 0; clip just for plotting.
    finite = np.isfinite(y)
    if np.any(finite):
        cap = np.nanquantile(y[finite], 0.995)
        y_plot = np.clip(y, 0, cap)
    else:
        y_plot = y

    fig.add_trace(go.Scatter(x=x, y=y_plot, mode="lines", name=f"c={c}"))

fig.update_layout(
    title="dweibull PDF shapes (clipped near the spike for c<1)",
    xaxis_title="x",
    yaxis_title="pdf(x)",
)
fig.show()

# CDF shapes for the same c values
x = np.linspace(-4, 4, 2000)

fig = go.Figure()
for c in c_values:
    fig.add_trace(go.Scatter(x=x, y=dweibull_cdf(x, c), mode="lines", name=f"c={c}"))

fig.update_layout(
    title="dweibull CDF shapes",
    xaxis_title="x",
    yaxis_title="cdf(x)",
)
fig.show()

6) Derivations#

(a) Expectation#

For the standardized distribution \(Z\) the PDF is symmetric: \(f(z)=f(-z)\). Provided \(\mathbb{E}[|Z|]<\infty\) (true for all \(c>0\)),

\[ \mathbb{E}[Z] = \int_{-\infty}^{\infty} z f(z)\,dz = 0. \]

With location, \(X=\text{loc}+\text{scale}Z\), we get \(\mathbb{E}[X]=\text{loc}\).

(b) Variance via the Gamma function#

Compute the absolute moment for \(k> -c\):

\[ \mathbb{E}[|Z|^k] = \int_{-\infty}^{\infty} |z|^k \frac{c}{2}|z|^{c-1}e^{-|z|^c}\,dz = c\int_0^\infty z^{k+c-1}e^{-z^c}\,dz. \]

Substitute \(u=z^c\Rightarrow z=u^{1/c}\) and \(dz=\tfrac1c u^{1/c-1}du\):

\[ \mathbb{E}[|Z|^k] = c\int_0^\infty u^{(k+c-1)/c} e^{-u}\,\frac1c u^{1/c-1}\,du = \int_0^\infty u^{k/c} e^{-u}\,du = \Gamma\!\left(1+\frac{k}{c}\right). \]

For the variance, take \(k=2\):

\[ \operatorname{Var}(Z)=\mathbb{E}[Z^2]=\Gamma\!\left(1+\frac{2}{c}\right). \]

(c) Likelihood and log-likelihood#

Assume i.i.d. observations \(x_1,\dots,x_n\) from the standardized model (\(\text{loc}=0\), \(\text{scale}=1\)). The likelihood for \(c\) is

\[ L(c) = \prod_{i=1}^n \frac{c}{2}|x_i|^{c-1}\exp(-|x_i|^c). \]

The log-likelihood is

\[ \ell(c)=n\ln\frac{c}{2} + (c-1)\sum_{i=1}^n \ln|x_i| - \sum_{i=1}^n |x_i|^c. \]

Differentiating gives a score equation (no closed-form MLE in general):

\[ \ell'(c)=\frac{n}{c} + \sum_{i=1}^n \ln|x_i| - \sum_{i=1}^n |x_i|^c\ln|x_i|. \]

This is typically solved numerically (as SciPy does under the hood).

# Example: visualize the log-likelihood over c (standardized case)

x_data = dweibull.rvs(0.8, size=4000, random_state=rng)


def loglike_c(c: float) -> float:
    if c <= 0:
        return -np.inf
    return float(np.sum(dweibull_logpdf(x_data, c)))


c_grid = np.linspace(0.2, 4.0, 200)
ll = np.array([loglike_c(c) for c in c_grid])

fig = px.line(x=c_grid, y=ll, labels={"x": "c", "y": "log-likelihood"}, title="Log-likelihood vs c")
fig.show()

c_hat = float(c_grid[np.argmax(ll)])
c_hat

0.7919597989949749

7) Sampling & Simulation (NumPy-only)#

Inverse transform sampling#

Using the PPF, we can sample with a single uniform random variable, but an even simpler implementation uses the sign + magnitude story:

  1. Sample \(U\sim\text{Unif}(0,1)\) and set \(Y = (-\ln U)^{1/c}\). (This is Weibull sampling.)

  2. Sample an independent sign \(S\in\{-1,+1\}\) with \(\mathbb{P}(S=1)=1/2\).

  3. Return \(X = \text{loc} + \text{scale}\cdot S Y\).

This is NumPy-only and avoids SciPy entirely.

def dweibull_rvs_numpy(
    c: float,
    loc: float = 0.0,
    scale: float = 1.0,
    size: int | tuple[int, ...] = 1,
    rng: np.random.Generator | None = None,
):
    # Draw random samples from dweibull using NumPy only.
    _validate_params(c, scale)
    if rng is None:
        rng = np.random.default_rng()

    u = rng.random(size)
    u = np.clip(u, np.finfo(float).tiny, 1.0)

    # Magnitude ~ Weibull(shape=c, scale=1): Y = (-log U)^(1/c)
    y = np.power(-np.log(u), 1.0 / c)

    # Random sign
    s = np.where(rng.random(size) < 0.5, -1.0, 1.0)

    return loc + scale * s * y


# Quick validation: sample moments vs theory
c_sim = 0.8
n = 100_000
samples = dweibull_rvs_numpy(c_sim, size=n, rng=rng)

mean_mc = float(np.mean(samples))
var_mc = float(np.var(samples))

mean_th, var_th, *_ = dweibull_theoretical_stats(c_sim)
(mean_mc, var_mc, float(mean_th), float(var_th))

(-0.014103101477437954, 3.3469432497508005, 0.0, 3.323350970447843)

8) Visualization#

We’ll visualize three things for a chosen parameter set:

  • the PDF

  • the CDF

  • Monte Carlo samples compared to the theoretical PDF

c_vis = 0.8
loc_vis = 0.0
scale_vis = 1.2

x = np.linspace(-5, 5, 3000)
pdf_np = dweibull_pdf(x, c_vis, loc=loc_vis, scale=scale_vis)
cdf_np = dweibull_cdf(x, c_vis, loc=loc_vis, scale=scale_vis)

pdf_sp = dweibull.pdf(x, c_vis, loc=loc_vis, scale=scale_vis)
cdf_sp = dweibull.cdf(x, c_vis, loc=loc_vis, scale=scale_vis)

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=pdf_np, mode="lines", name="NumPy pdf"))
fig.add_trace(go.Scatter(x=x, y=pdf_sp, mode="lines", name="SciPy pdf", line=dict(dash="dash")))
fig.update_layout(title="PDF: NumPy vs SciPy", xaxis_title="x", yaxis_title="pdf(x)")
fig.show()

fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=cdf_np, mode="lines", name="NumPy cdf"))
fig.add_trace(go.Scatter(x=x, y=cdf_sp, mode="lines", name="SciPy cdf", line=dict(dash="dash")))
fig.update_layout(title="CDF: NumPy vs SciPy", xaxis_title="x", yaxis_title="cdf(x)")
fig.show()

# Monte Carlo histogram vs theoretical PDF
n = 20_000
x_samp = dweibull_rvs_numpy(c_vis, loc=loc_vis, scale=scale_vis, size=n, rng=rng)

hist = px.histogram(
    x=x_samp,
    nbins=80,
    histnorm="probability density",
    opacity=0.5,
    title="Samples (histogram) vs theoretical PDF",
    labels={"x": "x"},
)

pdf_line = go.Scatter(x=x, y=pdf_sp, mode="lines", name="theoretical pdf")

fig = go.Figure(hist.data)
fig.add_trace(pdf_line)
fig.update_layout(xaxis_title="x", yaxis_title="density")
fig.show()

9) SciPy Integration (scipy.stats.dweibull)#

SciPy offers a full suite of distribution methods:

  • pdf, logpdf, cdf, ppf, rvs

  • stats for moments

  • entropy

  • fit for parameter estimation (numerical)

We’ll generate synthetic data, fit the parameters, and overlay the fitted PDF.

# Synthetic data
c_true, loc_true, scale_true = 0.9, -0.3, 1.4
x_obs = dweibull.rvs(c_true, loc=loc_true, scale=scale_true, size=5_000, random_state=rng)

# Fit all parameters
c_fit, loc_fit, scale_fit = dweibull.fit(x_obs)
(c_fit, loc_fit, scale_fit)

(0.8918390627792722, -0.2964891361570203, 1.3802126275005553)

x_grid = np.linspace(np.quantile(x_obs, 0.001), np.quantile(x_obs, 0.999), 2000)

pdf_true = dweibull.pdf(x_grid, c_true, loc=loc_true, scale=scale_true)
pdf_fit = dweibull.pdf(x_grid, c_fit, loc=loc_fit, scale=scale_fit)

hist = px.histogram(
    x=x_obs,
    nbins=80,
    histnorm="probability density",
    opacity=0.4,
    title="SciPy fit: true vs fitted PDF",
    labels={"x": "x"},
)

fig = go.Figure(hist.data)
fig.add_trace(go.Scatter(x=x_grid, y=pdf_true, mode="lines", name="true pdf"))
fig.add_trace(go.Scatter(x=x_grid, y=pdf_fit, mode="lines", name="fitted pdf", line=dict(dash="dash")))
fig.update_layout(xaxis_title="x", yaxis_title="density")
fig.show()

10) Statistical Use Cases#

(a) Hypothesis testing (goodness-of-fit)#

A common workflow is:

  1. fit parameters

  2. test whether the fitted distribution plausibly generated the data

A classic tool is the Kolmogorov–Smirnov (KS) test.

Caution: if you fit parameters on the same data you test, the KS p-value is only approximate (the null distribution changes). Still, it’s a useful diagnostic.

(b) Bayesian modeling#

The log-likelihood \(\ell(c)\) makes it easy to do Bayesian inference for \(c\) with a prior (e.g., Gamma prior). We’ll do a simple grid posterior example with known loc=0, scale=1.

(c) Generative modeling#

You can use dweibull as a drop-in noise distribution (especially for \(c\le 1\)) to generate data with heavier tails than a Gaussian and a different near-zero behavior.

# (a) KS test using fitted parameters (approximate when parameters are estimated)
D, p_value = kstest(x_obs, "dweibull", args=(c_fit, loc_fit, scale_fit))
(D, p_value)

(0.015158632683738738, 0.19873182866725814)

# (b) Simple Bayesian inference for c (assuming loc=0, scale=1 known)

# Generate standardized data
c_true_bayes = 0.75
x_bayes = dweibull.rvs(c_true_bayes, size=1500, random_state=rng)

c_grid = np.linspace(0.2, 4.0, 600)

# Gamma prior on c: shape α, rate β
alpha, beta = 2.0, 1.0
log_prior = (alpha - 1) * np.log(c_grid) - beta * c_grid  # constants omitted

log_like = np.array([np.sum(dweibull_logpdf(x_bayes, c)) for c in c_grid])
log_post = log_like + log_prior
log_post -= np.max(log_post)
post_unnorm = np.exp(log_post)

# Normalize
Z = np.trapz(post_unnorm, c_grid)
post = post_unnorm / Z

# Posterior mean and MAP
c_map = float(c_grid[np.argmax(post)])
c_mean = float(np.trapz(c_grid * post, c_grid))

# 95% credible interval via numerical CDF
cdf = np.cumsum((post[:-1] + post[1:]) / 2 * np.diff(c_grid))
cdf = np.concatenate([[0.0], cdf])

c_lo = float(np.interp(0.025, cdf, c_grid))
c_hi = float(np.interp(0.975, cdf, c_grid))

(c_true_bayes, c_map, c_mean, (c_lo, c_hi))

(0.75,
 0.7455759599332219,
 0.7472222579871374,
 (0.7182876018799911, 0.7762668748027535))

fig = go.Figure()
fig.add_trace(go.Scatter(x=c_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=c_true_bayes, line_dash="dash", line_color="black", annotation_text="true c")
fig.add_vline(x=c_map, line_dash="dot", line_color="royalblue", annotation_text="MAP")
fig.update_layout(title="Posterior over c (loc=0, scale=1 assumed)", xaxis_title="c", yaxis_title="density")
fig.show()

# (c) Generative modeling example: a smooth signal + different noise models

t = np.linspace(0, 1, 400)
y_true = np.sin(2 * np.pi * t)

sigma = 0.25
noise_gauss = rng.normal(0.0, sigma, size=t.size)
noise_dw = dweibull_rvs_numpy(c=0.7, scale=sigma, size=t.size, rng=rng)

y_gauss = y_true + noise_gauss
y_dw = y_true + noise_dw

fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=y_true, mode="lines", name="true signal", line=dict(color="black")))
fig.add_trace(go.Scatter(x=t, y=y_gauss, mode="markers", name="Gaussian noise", opacity=0.6))
fig.add_trace(go.Scatter(x=t, y=y_dw, mode="markers", name="dweibull noise (c=0.7)", opacity=0.6))
fig.update_layout(title="Same signal, different noise distributions", xaxis_title="t", yaxis_title="y")
fig.show()

# Compare residual distributions
residuals = {
    "Gaussian": noise_gauss,
    "dweibull (c=0.7)": noise_dw,
}

fig = px.histogram(
    x=np.concatenate(list(residuals.values())),
    color=np.repeat(list(residuals.keys()), repeats=[t.size, t.size]),
    nbins=70,
    barmode="overlay",
    histnorm="probability density",
    opacity=0.5,
    title="Residual distributions",
    labels={"x": "residual"},
)
fig.show()

11) Pitfalls#

  • Parameter validity: \(c>0\), scale>0. Invalid values should error early.

  • Zeros in data: for \(c>1\), the PDF at 0 is exactly 0; if your data contains many exact zeros (rounding/quantization), the likelihood can behave strangely. For \(c<1\), the PDF is infinite at 0, so exact zeros can dominate fits.

  • Bimodality for \(c>1\): this is often unexpected if you think of the model as “noise around 0”.

  • MGF nonexistence: for \(c\le 1\), the MGF does not exist for all \(t\) (Laplace has a finite strip; \(c<1\) diverges for any nonzero \(t\)).

  • Numerical stability: prefer logpdf in optimization; for large \(|x|^c\) the PDF underflows to 0 (fine), but products of PDFs can underflow without logs.

  • Fitting: dweibull.fit is numerical and can be sensitive; consider fixing loc if you know the center, or providing good initial guesses in custom optimization.

12) Summary#

  • dweibull is a continuous, symmetric distribution on \(\mathbb{R}\) with shape parameter \(c>0\).

  • Its PDF \(\propto |x|^{c-1}e^{-|x|^c}\) creates three regimes: spike at 0 (\(c<1\)), Laplace (\(c=1\)), and bimodal (\(c>1\)).

  • Key identity: \(\mathbb{E}[|Z|^k]=\Gamma(1+k/c)\), giving closed-form variance and kurtosis.

  • Sampling is simple via inverse transform / sign + Weibull magnitude, and SciPy provides pdf/cdf/rvs/fit utilities.

  • In practice, dweibull can be a flexible tool for diagnostics, Bayesian inference over shape, and generative noise modeling when Gaussian assumptions are not appropriate.

References

  • SciPy docs: scipy.stats.dweibull (notes include the defining PDF)

  • Standard Gamma function identities for Weibull moments

Contents